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We study the production of gravitational waves from cosmic domain walls created during phase 
transition in the early universe. We investigate the process of formation and evolution of domain 
walls by running three dimensional lattice simulations. If we introduce an approximate discrete 
symmetry, walls become metastable and finally disappear. This process might occur by a pressure 
difference between two vacua if a quantum tunneling is neglected. We calculate the spectrum of 
gravitational waves produced by collapsing metastable domain walls. Extrapolating the numerical 
results, we find that the signal of gravitational waves produced by domain walls whose energy scale 
is around 10^'^-lO^^GeV will be observable in the next generation gravitational wave interferometers. 
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I. INTRODUCTION 

The spontaneous symmetry breaking is one of the most important ingredients of the modern particle physics. Not 
only it describes the unification of forces acting on the elementary particles, but also it has rich consequences in 
the early universe. One is the occurrence of primordial phase transitions followed by the formation of topological 
defects [1]. Especially, if there was the spontaneous breaking of a discrete symmetry, the two dimensional surface-like 
defects, called domain walls, would be formed in the early era of the universe (for a review of domain walls and 
other topological defects, see [2 ). Although discrete symmetries naturally arise in many particle physics models, 
the existence of stable domain walls is disfavored by cosmological considerations. Indeed, if domain walls survived 
until the present time, the energy density of walls would eventually come to dominate the total energy density of the 
universe, causing the fast expansion of the universe to affect the galaxy formation or primordial nucleosynthesis [3]. 
Moreover, the presence of domain walls would cause the excessive anisotropy in the cosmic microwave background 
observed today, which rules out the existence of stable domain walls with the symmetry breaking scale r] >lMeV [4]. 

However, it is still possible to consider the existence of unstable domain walls which decay early enough not to lead 
cosmological disasters. One way to introduce the instability of walls is to make discrete symmetry only approximate 
and tilt the potential so as to lift the degeneracy of vacua. This might be raised by nonrenormalizable operators 
suppressed by the cutoff scale (~ Planck scale) [5 . The idea of an approximate discrete symmetry as a mechanism 
to remove the domain wall problem has been studied by several authors [3l[6]-[8]. Also, there is another possibility to 
create unstable domain walls due to a postinflationary nonthermal phase transition [9l [10] . 

In this paper we investigate the production of gravitational waves from the annihilating process of unstable domain 
walls. There are early works studying the bubble collisions in the early universe [11 and the possibility of generating 
gravitational waves p!2HT6] . In these works, gravitational waves are considered to be generated by the collision of 
spherically expanding bubbles, which are created during first order phase transitions (i.e. the quantum tunneling 
from the false vacuum to the true vacuum). Instead, in this paper we assume that the transition from one vacuum 
to another proceeds via a realignment of the classical field rather than a quantum tunneling. The production of 
gravitational waves from such a process was originally calculated in p/7j. However, in (TTJ the radiated energy is 
calculated only in the simple configuration such as situations with axial symmetry, and only the total intensity of 
gravitational waves was estimated. Instead, our purpose is to study the production of gravitational waves in more 
realistic situations (i.e. to follow from phase transition and to consider generic configuration of the scalar field), and 
calculate the spectrum of gravitational waves. 

The gravitational wave is the powerful tool for studying cosmology and high energy physics (for reviews, see p^H2Q] ). 
There are various ongoing and planned experiments on the gravitational waves. As for the ground-based experiment, 
LIGO [21] is running, and LCGT [22 is planned. And space-borne interferometers such as LISA [23 , BBO [M], and 
DECIGO [25 are planned to be launched in the future. Furthermore, it has been started to plan the ground-based 
interferometer with improved sensitivity, "Einstein Telescope (ET)" [26| in Europe recently. In this paper we discuss 
whether the stochastic gravitational wave backgrounds produced by domain walls are detectable in these experiments. 
If it is detectable and the predicted form of the spectrum is distinctive, we expect that the gravitational wave spectrum 
from domain walls can be another probe for high energy physics. One example is in the theory with supersymmetry. 
One may be able to measure the parameter of the theory such as the mass of the gravitino, fermionic partner of the 
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graviton, by gravitational waves from annihilation of Z2N domain walls [2T. 

The dynamics of domain walls is not easy to investigate because of the nonlinear nature of the evolution equations. 
Thus it is necessary to use numerical calculations to follow the evolution of domain walls with arbitrary structures. The 
first complete numerical study was performed by Press, Ryden, and Spergel [28 by using the modified field equation 
for the scalar field (called the PRS algorithm). Also, Kawano [29 simulated the evolution of domain walls by using 
the thin wall effective equation of motion (similar to the Nambu-Goto equation [30l |3T] for cosmic strings). Later, 
several authors investigated further by up-dated high resolution simulations [32"^3T and the numerical simulations 
for unstable domain walls are performed in [TOl [35]. However, although these numerical simulations successfully 
demonstrate the macroscopic evolution of domain wall networks, they put an unphysical assumption about the width 
of a domain wall, and we cannot use their computational method since this method may invalidates the estimation of 
the energy of gravitational waves (see Section IV A). Therefore, we solve the field equations directly rather than rely 
on the previous methodologies. This gives us severe technical restrictions in choosing parameters of the simulations. 

This paper is organized as follows. In Section |TI1 we shortly review the dynamics of domain wall networks. In 



Section III, we describe the formalism to calculate classical gravitational wave spectra radiated by the scalar field 



in the systematic way. Then, after a short discussion about the difficulty in numerical simulations, we present the 



results of our numerical simulations in Section IV We combine the results of numerical calculations with analytic 
estimations, and discuss the observability of the spectrum by extrapolating them in Section jV] Finally we conclude 
in Section [yll 



II. DOMAIN WALL DYNAMICS 



A. The Model of Domain Walls 



We consider the model of the real scalar field with the potential 

v{<p) = ^{4>' - v' f + ev4> I^l4>' - v') + It'^\ (1) 

The first term is the usual double well potential. In the absence of the second term, this model has a discrete Z2 
symmetry corresponding to the transformation (j) —cj). The second term is introduced in order to lift the degeneracy 
of the two vacua and make the discrete symmetry approximate. In this model, the scalar field is assumed to be thermal 
equilibrium at the first stage of their evolution, and we introduce the finite temperature correction in the last term 
of the potential. When finite temperature effects become negligible, the potential has two minima (j) = ±7^ and the 
formation of domain walls occurs. The difference of the energy density between two minima is A = Aerj^/S and the 
height of the potential barrier is AV ^ Xrj^/A. The dimensionless constant parameter e is called "bias", which makes 
domain walls unstable. 

In the numerical simulations which we will describe later, we solve the classical filed equation with the potential 
given by Eq. ([T]). This corresponds to the fact that we consider only a quiescent second order transition. If there 
are more violent first order transitions via a quantum tunneling, the gravitational wave signature might be enhanced 
significantly. Such a quantum treatment is, however, beyond the scope of this paper. 

In the following two subsections, we briefiy review the property of the evolution of domain wall networks. 



B. Scaling Solutions 



The well known property of the evolution of domain wall networks is that there exists the phase of "scaling", in 
which the typical length scales such as the wall curvature radius R and the distance of two neighboring walls L are 
given by Hubble radius, 

H-^ - t. (2) 

This property is originally studied in [28] numerically, and by other groups in I32ll35] . However, the analytic 
modeling of the scaling solution is found in [36^-^. 

When domain wall networks are in scaling regime, the energy density of domain walls is given by 

- aR^/R^ - cr/t, (3) 

where a = 2a/2A^^/3 is the tension of the wall. Eq. ^ is also expressed as 

A/V(xr-\ (4) 
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where AjV is the comoving area density of the wah, and r is the conformal time: dr = dija. If we assume radiation 
dominated universe, AjV is proportional to which is a criterion to distinguish whether walls are in scaling 

regime. 



C. The Role of the Bias 



There are two effects caused by the bias in the potential. One is the bias of two vacua chosen during the phase 
transition, and another is the volume pressure which shrinks the false vacuum domains and eliminates domain walls. 
The cosmological scenario with biased discrete symmetry is investigated in [3 and here we follow the analysis of it. 

Let us denote the probability of having a scalar field fiuctuation during the phase transition end up in plus vacuum 
as p+ and in minus vacuum as p_ . If there exists the bias, the ratio of these two probabilities is given by 

exp — — . (5) 



p- '^y T 

Here, AF = A x ^^^^^ is the difference of the free energy between two vacua, and ^corr is the correlation length of the 
scalar fields. If we estimate T as Ginzburg temperature [Tl'3|, given by Tq — AV x ^^^^^ ~ Xrj'^^^^^^/A, the above 
equation can be written as 

^=exp(-^)^exp(-f I). (6) 

Therefore, there exists the bias between two vacua if e 7^ 0. The spatial distribution of two vacua after phase transition 
has been studied by using percolation theory [39|. It has been shown that if the probability p+ (or p_) is bigger than 
a critical value an infinite plus (minus) cluster appears in the space. If we approximate the space as a three 
dimensional cubic lattice, Pc = 0.311 [39 , and by requiring both p+ and p- are bigger than pc^ we obtain 

e < 0.15A. (7) 

If the condition given by Eq. ^ is satisfied, the infinite size of domain walls would appear after phase transition. 

After the formation of domain walls, there are two forces acting on the walls. The first force is a surface tension, 
which is given hy pt — cr/R — 2V2Xr]^ /3R (as a pressure). It smooths out the small scale structure on the wall and 
straightens the wall up to the horizon scale. The second force is a volume pressure, which is given 

It accelerates the wall against the false vacuum regions due to the energy difference between the two vacua and collapses 
the wall. The evolution of domain wall networks depends on when and which of the two forces dominate. The volume 
pressure dominates when pv — Pr is satisfied, and the corresponding wall curvature radius is 

Rc^^iev)-'. (8) 

If we require that the wall has been straighten up to the horizon scale before the volume pressure becomes effective, 
we obtain the condition R > H~^(tc) where R is given by Eq. ^ and H(tc) is the Hubble parameter at the phase 
transition (T = Tc = 2r]). This condition gives an upper bound for e: 

e<4.7x(A5,)V2_|_, (9) 

where is the number of relativistic degrees of freedom at Tc and Mp is the Planck mass. On the other hand, a 
requirement that domain walls do not dominate the energy density of the universe gives a lower bound for e. When 
domain wall networks are in scaling regime, the energy density of domain walls is estimated as py^ cf/H~^ :^ 
\/2\rf /?>t. If we assume the radiation dominated universe, the total energy density of the universe is given by 
p 3Mp/327rt^ . The domain walls start to dominate the universe {pw/p — 1) when t = twDi where 

Meanwhile, the typical time scale of the collapse of the wall network is given by the condition H~^{tdec) = 2tdec — R 
where R is given by Eq. ([8|: 

idec^^y|(e^)-^ (11) 
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Therefore, the condition that domain wahs collapse before wall domination is tdec < twD^ from which we obtain 

167r / T] Y 



(12) 



To sum up, the domains with infinite size appear if the condition given by Eq. ^ is satisfied, and the domain walls 
collapse before wall domination if e satisfies Eqs. ^ and (12). The decay time of the domain walls is given by Eq. 



III. CALCULATION OF GRAVITATIONAL WAVES 



In this section, we describe how to calculate the spectrum of the stochastic gravitational wave backgrounds produced 
by domain walls. We use the formalism of "Green's function method" developed by Dufaux et al. [40l[4T]. Originally 
this formalism was applied to the production of gravitational waves from preheating [42", after infiation. However, 
this framework is applicable to general scalar fields in expanding universe, especially the Higgs field forming domain 
walls considered here. In the following, we summarize the results of [40]. 

Assuming a spatially fiat Friedmann-Robertson- Walker background, gravitational waves are represented by the 
spatial metric perturbation 

ds^ = a^{T)[-dr^ + {Sij + hij)dx'dx^], (13) 
where hij satisfies the condition dihij = h] = 0. The evolution of hij is described by the linearized Einstein equation 

hij + 3Hh,j - -^h,j = 16^GT^^, (14) 

where is the transverse-traceless part of the stress-energy tensor, which is computed by applying the projection 
operator in the momentum space 

(r,k) = A,,-fc,(^)T,,(r,k) = A,,-fc,(^){Sfc(/)Sz(/)}(r,k), (15) 

KjMih = P^k{k)PJl{k)-^P,J{k)Pkl{k). (16) 
P{k) = Sij-hkj, (17) 

where k = k/|k| and {9/c^9^0}(r, k) is the Fourier transform of 9/c^(r, x)9^^(r, x). 

Suppose that the source term is nonzero during the interval < r < r/. In this case, the equation of motion (14) 
is solved by using the Green's function in < r < r/, and the solution can be matched to the source- free solution in 
r > Tf. In the radiation dominated universe, the result is [40] 

hij{r,k) = A^_^-(k)sin[^(r - r/)] + Bij{k) cos[k{r - Tf)] for r > r/, (18) 

where 



Aj{k) = — I dr'cos[k{rf-r')Hr')T^^^{r',k), 
i?,,(k) = ^ dr' sin[fc(r/ - r')]a(r')7;7(r', k), 



(19) 



and hij = a{r)hij. 

The energy density of gravitational waves is given by the spatial average of the product of hij (see e.g. [T9]) 

1- • 1 1 f d^k - 

= ^('^«.(^.-)^«.(^.-)) = ^^.v J J2^'^^^''^'^^''^^ (20) 

where a prime denotes a derivative with respect to the conformal time r and V is the comoving volume. Strictly 
speaking, there are other contributions for the right hand side of Eq. (20): one is cross terms between hij and h[j 
which turn out to vanish under the reality condition of hij{r^x.)^ and another is the term proportional to which 
is ignorable if we consider only the modes inside the horizon at present time k/ao ^ Hq. Substituting Eq. (18) into 
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Eq. (20) and taking the time average over a period of the oscihations caused by the sine and cosine terms in Eq. (18) 
to ehminate these oscihations, we find 

47rG 



a^{t)V 



(21) 



The present spectrum is represented by the fraction of energy density against the total energy density of the universe 
per logarithmic frequency interval 



Pc,o d\nf 



Pc,0 



gwV 

d\nk 



(22) 



where / is the frequency of gravitational waves, Pc,o is the critical energy density today, and h is the renormalized 
Hubble parameter: Hq = 100/ikmsec~^Mpc~^. Taking account of the fact that energy density of gravitational waves 
redshifts as /5gw '-^ ^ 

, and from Eqs. rt21k and (^22L we obtain 



1/3 



PR{ti) 2^ 



where 



y" cZfifc ^ I ^ ^ dr' cos(fcT')a(T')T^^(T', k) + J^^ dr' sm{kT')a{T')T^/ {t' , k) | . 



(23) 



(24) 



Here go = 3.36 is the number of radiation degrees of freedom today, ^nh'^ =4.15 x 10~^ is the ratio of the radiation 
energy density to the critical energy density today, and pR{ti) is the radiation energy density at the initial time. 
is a unit vector representing the direction of k and dQk = d cos Odcj). We set the scale factor at the initial time as 
a{ti) = 1. Meanwhile, the frequency observed today is given by the comoving momentum redshifted by a{to)~^ 



f = 



go 



27ra(to) 27r y^r* 



1/3 



To 



(25) 



where Tq = 2.725K is the temperature of the universe observed today and Ti is the initial temperature (to be specified 
later). 

The computational strategy is as follows: first, we obtain the time evolution of scalar fields in the real space (j){t^ x) 
by the lattice simulation, then we compute the TT projected stress-energy tensor as Eqs. (15)-(17), and finally we 
increment the time integral in Eq. (24) to obtain the spectrum by using Eqs. (23) and (25). 



IV. NUMERICAL SIMULATION 



A. Problem in Numerical Study and Choice of Parameters 



The evolution equation of domain walls is highly nonlinear and we are forced to use numerical analysis. However, 
there is a technical problem in computing the evolution of domain walls numerically: one must consider two extremely 
different length scales simultaneously. One is the width of the wall ~ 77"^, which must be bigger than a lattice 
spacing throughout the numerical simulation in order to resolve domain walls by individual unit of lattice. But if 
we perform simulations in the comoving box, the physical lattice spacing expands as ex a(t), and becomes unable 
to resolve the wall in the finite simulation time scale. The other length scale is the Hubble horizon, which must be 
smaller than the simulation box in order to avoid unphysical effect caused by the finite size of the simulation box. 
Unfortunately, there is a large gap between these two length scales by a factor of the ratio between the mass scale of 
the scalar field and the Planck scale as 

S^-V-'« [^y-'-H-\ (26) 

Therefore, we can not simulate the physical evolution of domain walls for a long time for energy scale 77 sufficiently 
smaller than the Planck scale. 
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The standard PRS algorithm [28^ avoids this difficulty described above by keeping the comoving wall width constant 
in time, to maintain the width resolvable throughout a simulation. However, in this work we can not apply this method, 
since the energy of gravitational waves radiated by domain walls is estimated incorrectly due to the unphysical nature 
of the wall width. 

Hence we directly solve the physical equation of motion 



dV 



0, 



(27) 



with the potential given by Eq. ([T]). This restricts the time scale of the simulation and the model parameters. We 
must choose 77 close to the Planck scale due to the fact represented by Eq. (26), but if we choose 77 close enough to 



Mp, the onset of scaling is delayed in the unit of Hubble time (we choose time coordinate normalized by the initial 
Hubble time ti ^ H~^{ti)^ see Appendix Al) due to the degeneracy between the inverse mass scale r]~^ and the 
Hubble radius H~^. As the marginal possible parameters, we choose A = 1.0 and r] = 1.05 x lO^^GeV. We assume the 
radiation dominated era and the relativistic degree of freedom is constant given by = 100 during the simulation. 
The final time is set to be = 151ti where ti is initial time chosen so that the initial temperature is twice of the 
critical temperature of the phase transition {Ti = 2Tc = 4?], ti = tc/4:). We performed lattice simulations with 256^ 
points in the cases with 6 = 25 and b = 50, where b is the (comoving) box size in the unit of ti. In these parameters, 
the resolution of the wall width is well maintained and the horizon scale is smaller than the simulation box even at 



the final time (see Appendix A 1 ). 

For biased domain walls, we choose the parameter e such that the conditions Eqs. ^ and (12) are satisfied: 



167r 



V 
Mp 



< e < 4.7 X (A^*) 



1/2 _n_ 
Mp 



(28) 



In the parameters described above, Eq. (28) corresponds to the condition 4.1 x 10 < e < 0.40. We run the 



simulations varying e in the range e = 0.015, 0.02, 0.025, 0.03, 0.035, and 0.04 which satisfy the condition Eq. (28). 
These values also satisfy the condition for an appearance of the infinite size of domains, Eq. (7|). 

We emphasize that the dynamical range of the simulations is rather small compared with realistic cosmological 
dynamical ranges. The actual dynamical range is 12 in conformal time. For the domain wall evolution, this is reduced 
to less than (151ti/30ti)^/^ ^ 2.2 since the result with e = do not approach scaling until t = 30ti (see next subsection 
and Fig. [2|. On the other hand, the typical time scale for the wall domination is 64 [according to Eq. ([lO|]. It is 
extremely difficult to improve dynamical ranges without using PRS algorithm. Future simulations with much larger 
dynamical ranges should confirm and improve the results of our current simulations. 



B. Evolution of Domain Wall Networks 

The typical evolution of stable (e = 0) and biased (e 7^ 0) domain wall networks is shown in Fig. [l] The initial 
thermal fiuctuations of the scalar field are damped around t = 30^^, and domain walls extended over the horizon 
scale are formed. In the case e = 0, the wall networks evolve with annihilation and reconnection to maintain scaling 
property. If we turn on the bias, instead, the region occupied by false vacuum after the phase transition is smaller 
than the unbiased case, and walls gradually collapse due to the volume pressure. 

To see more quantitatively, we calculate time evolution of the comoving area density of domain walls as shown 
in Fig. [2] (see Appendix A3 for a detail of the calculation). From this figure we can see that the scaling regime in 



which the area density scales as begins around t = 30ti for the case e = 0. In the 6 = 25 case, the scaling law 

breaks down around the late stage of the evolution. This may be caused by an unphysical effect of the finiteness of 
the simulation box, that is, the horizon scale become comparable to the size of the simulation box around the final 
time (the similar deviation from the scaling solution can be found in [28l [33^ 134 ] ). 

In the primary thermal stage t < 30t^, the correlation length of the scalar field is determined by mass scale of the 
scalar field ^corr ^ 'rn~^ ^77"^ rather than the Hubble horizon scale. In this stage, the value of the scalar field changes 
randomly in the short length scale given by ^corr, thus A/V remains nearly constant. Note that if bti ^ ^corr ^ ^ 1 
by increasing box size b with the fixed number of grids one can count more points where the scalar field change their 
sign, leading to increase the value of A/V. In fact, the value of A/V at the primary stage t < 30ti in the case 6 = 50 
is larger than that in the case 6 = 25, as shown in Fig. [2j 

As we turn on the bias, the value of A/V falls by a factor of 0(10-100) in the simulation time scale. The typical 



time scale of the decay of the area density is well described by the estimation Eq. (11 ) except a factor of 0{1). There 
are oscillations of A/V after the decay of domain walls. This can be interpreted as the radiation of the energy stored 
in walls as the scalar field oscillating around the true vacuum [35]. The difference in the amplitude of oscillations 





FIG. 1: The evolution of stable (e = 0: left) and unstable (e ■ 
region where the value of crosses the zero. 



0.02: right) domain walls. The white surface corresponds to the 



between the case with 6 = 25 and 6 = 50 may be due to the poor resolution in the simulation with h = 50. If the 
wavelength of the scalar field radiations produced by the decay of domain walls is comparable or smaller than the 
small scale cutoff (~ ^~^), the typical length scale over which the scalar field varies becomes shorter than the lattice 
spacing. This makes the value of A/V inaccurate because we estimate the value A/V as a sum of the link of points 
where (j) has different signs (see Appendix A3). In the case with 6 = 50, the lattice spacing dx is larger than that 
with 6 = 25 and it is likely to lose many links over which the scalar field changes its signs. Therefore, the oscillation 
of the A/V in the case with 6 = 50 is less violent than that with 6 = 25, as we see in Fig. |2j 

C. Spectrum of Gravitational Waves 



By using the method described in Section [Illj we compute the spectrum of gravitational waves produced by domain 
walls. The results are shown in Fig. [3j To take small b with fixed number of lattice points corresponds to increasing 
the spatial resolution. Therefore, the frequency range computed in the case with 6 = 25 is higher than that in the case 
with b = 50. The two figures shown in Fig. [3] seem to be different in spite of the fact that the two results are obtained 
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FIG. 2: The time evolution of the comoving area density of domain walls for various values of e. The left panel shows the result 
in the case with the box size b — 25 and the right panel shows that in the case with 6 = 50. 

from calculations with completely the same parameters except b. For example, there is a peak in high frequency edge 
around 2 x lO^^Hz of the spectrum with b = 50, but the corresponding peak is not found in the spectrum with b = 25. 
Hence this high frequency peak in 6 = 50 is thought to be unphysical noise due to the poor resolution for a small 
scale. Another difference is that the low frequency tail of the spectrum in 6 = 25 is steeper than that in b = 50. This 
is due to the finite volume effect which promotes walls to collapse around the final time in 6 = 25, as shown in Fig. 
[2] Because of the limitation in our computer memory, we can not perform the simulation with maintaining physical 
resolution over full range of gravitational wave spectra. From now on, therefore, we use the result of the case with 
6 = 50 and ignore the high frequency peak around 2 x lO^^Hz. 
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FIG. 3: The spectrum of gravitational waves from domain walls for various values of e. The left panel shows the result in the 
case with the box size 6 = 25 and the right panel shows that in the case with 6 = 50. 



In Fig. [4) we show the time evolution of the gravitational wave spectra for several values of e in the case b — 50. 
There are also the unphysical high frequency peak in those spectra, but we ignore them for the reason described 
above. By comparing Fig. [2] and Fig. |4]we can see that the spectrum of gravitational waves evolves correspondingly 
to the evolution of the scalar field. In the primary stage {t < 30^^), the thermally distributed scalar field produces 
the spectrum which has a peak around the frequency / ~ 10^^ Hz, and after the formation of scaling wall networks 
(t > bOti) there is an increase in gravitational wave amplitudes of low frequency modes. These low frequency modes 
turn out to correspond the horizon scale at the time when gravitational waves are produced. We will investigate it 
in detail in Section |Vl In the late stage, if e = 0, the amplitudes of the gravitational waves keep growing in both 
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high frequency of the smah scale and low frequency of the horizon scale. On the other hand, if e 7^ 0, the growth in 
low frequency terminates when walls begin to collapse, and only high frequency modes grow in the late stage. The 
resulting spectrum has the low frequency edge and the high frequency peak, and gradually increases in intermediate 
modes. 
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FIG. 4: The time evolution of the gravitational wave spectra for various values of e (in the case b = 50). The spectra are shown 
from the time t = llti (pink) to t = 151ti (green) with the interval At = 20ti. 



V. ANALYTIC ESTIMATION AND PERSPECTIVE FOR OBSERVATIONS 



In this section we discuss the physical interpretation of the numerical results. We check whether the results of 
the numerical simulation match the naive analytic estimations. Finally we extrapolate the results into the general 
parameter space and estimate the amplitude and the frequency of gravitational waves observed today for various 
parameters. In particular we discuss whether the spectrum is observable in the future gravitational wave experiments. 



A. Physical Interpretation of the Spectrum 



There are two mechanisms relevant for production of gravitational waves from domain wall dynamics. One is due 
to the relativistic motion of walls in the universe, and another is due to the interactions of domain wall networks 
such as collision, separation, and reconnection of walls. Considering the scaling property of domain walls, the typical 
scale relevant for the former is the horizon scale at the time when gravitational waves are produced. Meanwhile, the 
typical scale relevant for the latter is the small scale of the structure on the domain wall. As long as domain walls are 
in the scaling regime, both of these two mechanisms are efficient, and the resulting spectrum extends over a broad 
frequency band between the low frequency of the horizon scale and the high frequency of the small scale structure of 
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the wall, as shown in Fig. [3] Instead, the growth of the amplitude in low frequency modes terminates when domain 
walls begin to collapse if there is a bias. 

From the above considerations, we conjecture that the form of the spectra is determined as follows. There is a 
peak in the high frequency corresponding to the small scale structure of the wall (i.e. the wall width). The spectrum 
gradually decreases toward the low frequency, and there is an 'edge' in the low frequency corresponding to the horizon 
scale. For causality reasons, there are no modes beyond the horizon scale. 

In the following, we estimate the amplitude and the frequency of gravitational waves analytically in order to justify 
the above considerations. 



B. Analytic Estimation 



First we estimate the peak frequency. The peak is given by the interaction of domain wall networks and we naively 
estimate that the peak frequency corresponds to the width of the domain wall. One can think of this process as 
the collision of wave packets of classical scalar fields whose wavelength is about Apeak ^ ^ (A^/^?])"^. Then, the 
physical wave number of the peak of gravitational waves at the time when they are produced is k^eo,]^/ a{t^) = 
27r/Apeak = 27rA^/^7^, where /Cpeak is the comoving wave number corresponding to the peak frequency. Taking redshift 
into account, the peak frequency observed today is obtained as 



^peak _ ( CL{t^) ^ ^1/2 



Setting to the final time of the simulation, 151t^, the estimation of the peak frequency gives 3.5 x lO^^Hz, which 
agrees with our numerical results except the factor of 0{1). Note that if we assume the entropy conservation, the 
dependence /Cpeak cx is canceled by the factor a{ti)/a{to) ex 1/r]. Therefore, generically the peak frequency is located 
in the high frequency of order lO^^Hz regardless of the value of rj. 

Next, we estimate the amplitude of gravitational waves observed today. The energy of gravitational waves produced 
by domain walls is estimated as Eg^ ~ GM^y^/ R^^ where is a typical length scale over which the energy distributes 
and Mdw is the mass scale of the domain wall. We do not specify the value of i^* since the result is not depend 
on R^ as shown shortly. We estimate M^w ~ o-Rl = V2Xr]^Rl/3^ and the energy density of gravitational waves 
at the production time is Pgw(^*) ^ ^gw/^* ~ {SX/9){r]/Mp)'^r]'^. By redshifting it, we find the amplitude of the 
gravitational waves observed today as 

[a{to)) Pc,o ^ ^ ^0 U J Vl-05 X lOi^Gevj \l^iu) ' 

which agrees with the peak amplitude ^^^^K^ ~ O(10~^) of the numerical results. 

We also estimate the location of the low frequency edge in the spectrum which corresponds to the horizon scale. 
The physical wave number corresponding to the horizon scale at time t is given by khor/ci{t) = 27TH{t) = ir/t. Hence 
the frequency observed today is given by 

/^^^ = = ( ^] 1 = 6.58 X 10^^ X (t/t,)-'/' ( \v7^) Hz. (31) 

27ra(to) \a{to) ) 2t ^ ' \im x IQi^GeV ) ^ ^ 

Notice that ti ex 7^~^ and hence /hor is independent of both r] and ti as it should be. For example, the low frequency 
edge at the final time of the simulation t = 151^^ becomes /hor — 5 x lO^Hz, which reproduces the location of the 
edge in the spectrum obtained by numerical simulations. 

Before moving to the next subsection, we discuss the origin of the spectrum produced during the primary stage 
{t < 30ti) of the scalar field evolution. In this stage, the scalar field fiuctuates randomly with correlation length 
roughly given by inverse of the temperature ^corr T~^, and gravitational waves can be produced via the quadratic 
interaction of the scalar field. We denote the amplitude of the scalar field fiuctuations by 5(j). From the equation of 
motion for the metric perturbation, Eq. ( p!4| ), a typical amplitude of gravitational wave hg with the comoving wave 
number k is estimated as k'^hg ~ 167rV^V^/Mp ~ 167r/c^(50^/M|,, thus we obtain hg ^ 50 x dcjp' /Mp. On the other 
hand, the ratio between the energy density of gravitational waves and the total energy density of the universe at the 
production time is estimated as (l^gw)* ~ Pgw/pnit^) ~ h'^/GpR{t^) ~ h'^/H'^{t^) ~ tlh'^ ^ /i^, where pR{t^) is the 
energy density of the radiation at t = and we assume the radiation dominated universe. In our simulations, the scalar 
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field initially fluctuates with amplitude dcj) Ti r] lO^^GeV, which gives (llgw)* ^ 10^ x {5(l)/MpY ^ 10"^ 
The amplitude of these gravitational waves observed today is flG] 

/100\ 

^^^h^ = 1.67 X 10-^ l-^j (l^gw)* - 10-^^ (32) 

This agrees with our numerical results (see the spectra of t = llt^ and 31^^ in Fig. [4|. The location of the peak 
(i.e. / ^ lO^^Hz) is almost unchanged from the early time to the flnal time. This is because the typical frequency 
of gravitational waves produced in the early time is roughly given by correlation length of the scalar fleld which is 
determined by the mass scale m ~ a/3A^, and this scale is nearly equal to the typical scale of the small structure on 
the domain wall ~ 



C. Fitting the Numerical Results 

The numerical results reveal that the gravitational wave spectra appear in rather high frequency band. This is 
partially due to the fact that we assume very high energy scale, r] lO^^GeV. For the case with lower 77, it is 
impossible to perform the numerical computation correctly because of the technical restriction described in Section 
|IV A| Instead, we assume the parameter dependence of analytic estimations obtained above, and flt the numerical 
factor into the results of numerical simulations to derive the formulae for the frequency and the amplitude of peak 
and edge in the gravitational wave spectra. In deriving these formulae, we check the e dependence of the numerical 



results. The difference in the bias e affects the time scale of the collapse of the domain wall networks [see Eq. (11)]. 
In the following we assume that the typical time of the production of gravitational waves is identical to the decay 
time of the domain walls tdec, and estimate e dependence of the spectrum. Then we see if our numerical results follow 
the dependence. 



Let us consider the peak frequency. Assuming the parameter dependence of Eq. ( 29 ) and replacing the numerical 
coefficient of 0{1) by that of the numerical result, we obtain 

1/3 / , X 1/2 



1/12 , ^ ,1/2 

Hz. (33) 



In the second line, we substitute t^ec given by Eq. ( pT] ) into t^. The difference in dependence between the flrst 
line and the second line comes from ti = (45Mp/167r^^^T/)i^^ ex g^^^'^r]~'^ which is determined by the Friedmann 



equation at the initial time. From Eq. (33), we expect that the location of the peak shifts to lower frequency for a 
larger value of e. This is because if e is large, domain walls collapse in the early stage of their evolution, and the 
gravitational waves redshifts for a longer time than the case of a small e. In fact, the peak in our numerical simulations 
shifts to the lower frequency for a larger value of e (see Fig. [3|. But if e is sufficiently large, the peak mixes with 
the primary spectrum produced at t < 30^^, therefore it is difficult to check precisely the e dependence described by 



Eq. (33). 



Next we consider the peak amplitude of the gravitational waves. Substituting o:^ tdec into Eq. (30), we obtain the 

parameter dependence of the amplitude l^gw/?'^ cx g^^^^ Xr]'^ (tdec/ti)'^ ^ g^^^^ \^ri'^e~'^ . So the amplitude is proportional 
to 6~^. In Fig. (5) we plot the e dependence of the peak amplitudes obtained by numerical results for h = 50, and we 
found that the property Vtg^K^ cx is indeed observed in our numerical results. By fitting the numerical coefficients 
into the line in Fig. [Sj we obtain the peak amplitude 

Note that the result with e = 0.015 slightly deviates from the line of in Fig. [sj There are some possibilities 
to explain this deviation. One explanation is that the deviation is due to the lack of the simulation time. If e is 
sufficiently small, the domain walls have not completely collapsed in the simulation time. Therefore, there would be 
growth of Vtg^h^ if we run the simulation for a longer time. Another possibility is that the relation Qg^h'^ cx is 
not exactly satisfied for small e. In fact, if we take the limit e ^ 0, Vtg^h^ diverges. Moreover, we must be careful 
about the contribution of the energy density of domain walls to the total energy density of the universe for the case 
with small e, which is ignored here. Unfortunately, we can not test these possibilities in our current capability of 
computers. 
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Q. 




FIG. 5: The relation between {^gwh^)peak and e obtained from the numerical results. The dotted line represents the fitting 
formula given in Eq. ( 34 ) . 



Now, we consider the intermediate frequency range between the high frequency peak and the low frequency edge. 
We parameterize the form of the spectrum in this range by Q.g^h? (x f^. In Fig. [6j we plot the value of n from 
numerical results in the case with the box size 6 = 50 for e =0.015, 0.02, 0.025, and 0.03. We also include the data 
with e = for the power law fitting in Fig. [6j The result of e = is not the spectrum which we would observe because 
the walls do not annihilate in the simulation time. However, we expect that the form of the spectrum with e <C 1 is 
well approximated by the result with e = if we assume that the low frequency modes of the gravitational waves are 
generated by the scaling evolution of the wall networks. Note that we omit the data for e = 0.035 and 0.04 since the 
walls collapse without following scaling phase sufficiently and the form of the spectrum significantly different from 
that of the small value of e. By fitting the numerical results, we obtain the frequency dependence as 

^gw/^' ^ /20e+0.08 fo^ f^^^^ f^^^^^ (35) 

For a large value of e, the domain walls immediately collapse after the formation, and there is little growth of the 
amplitude in the lower frequency modes corresponding to the horizon scale, resulting in a steeper spectrum. On the 
other hand, for a small value of e, domain walls are straightened up to horizon scale before they collapse. In this 
case, the difference in the value of Vt^-^h? between the peak and the edge is determined by the duration time of the 
collapse, which might be independent of e. Therefore, the value of the spectral index is independent of bias for e <C 1 



and then n 0.08. We emphasize that there might be large uncertainties in the parameterization in Eq. (35). If we 
write the frequency dependence as ftg^h'^ oc f^^+^ ^ the standard errors for the coefficients a and {3 are about 10% 
and 50%, respectively. These errors arise from the fact that there are only a few data points because of the lack of the 
dynamical range of the simulations. It will cause large uncertainties in the magnitude of gravitational waves relevant 
to the observations when we extrapolate the result in Section [VD| 

Before closing this subsection, we consider the edge in low frequency. We can not read the precise location of the 
edge from the numerical results because of the large statistical variance in the low frequency modes. Therefore, here 
we just derive the parameter dependence of the frequency and the amplitude for the low frequency edge analytically, 
instead of fitting numerical coefficients into the results of the numerical simulations as the previous discussions. 



Substituting t = t^ec into Eq. (31), we find 

/edge = /hor(tdec) = 5 X lO^" X £1/2^-1/4 1^^^ (_^)'^'hz. (36) 

For small e, the life time of domain walls becomes long, and the horizon scale at tdec becomes large. Then /edge shifts 
to the lower frequency, as shown in 6 = 50 case of Fig. [3j 



Also, by using Eq. (35), we can estimate the amplitude of gravitational waves at the edge as 



(^gw^ )edge — I ^ 1 (^gw^ )peak' l^'^) 

\ /peak / 
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FIG. 6: The relation between the spectral index n and the bias e. The dotted line represents the fitting formula given in 
Eq. (|35l. 



Substituting Eqs. (33), (34), and (36) into Eq. (37), we obtain 



-2 A -2 



((^gw/i^)edge ^ 3 X 10"-^'' X e-^X 



3g^^20e+0.08 

V 9^ 



v y 



(38) 



D. Forecasts for Observations 



In the previous subsection, we found the fitting formulae given by Eqs. (33)-(36) and (38) by combining the results of 



the numerical calculations and the analytic estimations. Now we use them to predict the gravitational wave spectrum 
for generic values of parameters. For simplicity we take A = 1.0 and = 100 in the following, then the spectrum 
depends on two parameters, the mass scale of the scalar field 77 and the bias e. Note that e must satisfy the condition 
given by Eq. ([28| so that our results are applicable. 

At a glance of the formulae (33)-([38|, we can say that r] determines the frequency of the spectrum and the strength 
of the amplitude, while e determines the amplitude and the bandwidth of the gravitational wave spectrum. The large 
value of T] results in the appearance of the gravitational wave spectrum in high frequency range and large flg^h'^. The 
small value of e results in the broad band of spectrum and also large Vt^^h^ . In Table |lj we show some examples of 
the value of various quantities estimated by using Eqs. (33)-(38). 



From Table [l| we see that if the value of r] is around lO'^^^-lO'^^GeV, the edge would appear in the frequency band 
observable in future gravitational wave astronomy projects such as ET, BBO, and DECIGO. On the other hand, the 
peak of spectra appears in too high frequency to be observable. In Fig. [7| we show examples of gravitational wave 
spectra from domain walls estimated by using the formulae (33)- (38). We note that the maximum value of Vtg^h^ at 
the peak is smaller than the bound which comes from Big Bang Nucleosynthesis (BBN) [TSl [48 . We also show the 
contours of the frequency and the amplitude of the edge in the parameter space of r] and e in Fig. [8) The parameter e 
is constrained by the condition given by Eq. (28), and the white region in Fig. [s] is allowed by this constraint. The red 
region denoted by "Not scaling" is the parameter space in which walls disappear before they straighten up to horizon 
scale. The gray region denoted by "Wall domination" is the parameter space in which the energy density of domain 
walls dominates the total energy density of the universe. Note that the red region and the gray region do not imply 
that parameters in these regions are ruled out. In fact, domain walls can decay before reaching the scaling regime 
if e is sufficiently large. Also, it might be possible that domain walls dominate the energy density of the universe if 
they disappear before the BBN epoch. However, we can use the formulae (33)-(38) only in the white region because 



in our analysis we assumed that domain walls straighten up to the horizon scale and disappear before they dominates 
the energy density of the universe. Especially, our numerical simulations are performed in the radiation dominated 
background, and we do not take into account the contribution of the energy density of domain walls. From Fig. [8| 
we can say that if r] is around lO^^-lO^^GeV, we would be able to measure the value of the scale of 77 and the bias 
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FIG. 7: The schematics of gravitational wave spectra from collapsing domain walls. We show three spectra for the case with 
(7y,e) = (10^2GeV,10-^^) (blue), (ry, e) = (10^°GeV, 10"^^) (red), and {r],e) = (10^°GeV, 10"^^) (green). We also roughly 
plot the sensitivity of planned detectors. Advanced LIGO 44 , LISA [23], ET [l6], BBO [24], and Ultimate DECIGO [45]. 
The dotted lines represent the constraints obtained from the observation of millisec pulsars [46l [47] (black) and Big Bang 
Nucleosynthesis [18. .48^ (green). 



e from the frequency and the amplitude of the edge of the gravitational wave spectra in future gravitational wave 
experiments. 

Before going to the conclusion, we comment on the accuracy of the estimation we made. There are large uncertainties 
in the determination of the frequency dependence in Eq. (35 ) because of the limited dynamical range of the simulations. 
As we noted before, the uncertainty of the spectral index n is about 50%. This uncertainty might become larger when 
we extrapolate the result for other values of e. Even if we neglect this increase in the uncertainties and only take the 
uncertainty in the numerical result into account, we expect that the uncertainty in the magnitude of (^gw^^)edge for 



TABLE I: The estimation of the amp litude and the frequency of the peak and the edge for various value of 77. We also show 
the condition for e given by Eq. ( [28| . The upper line in each row corresponds to the value in the case where e is taken as 
its minimum value, and the lower line corresponds to the value in the case where e is taken as its maximum value. Note that 
domain walls with a low energy scale such as 77 = lOOGeV and small e might be ruled out by the pulsar timing experiments 
(see Fig.|7|. 



V 


fpeah 


(l^gw^' )peak 


/edge 


(rigw^' )edge 


Condition for 


e 


lO'^GeV 


4.5 X lO'^Hz 
1.4 X 10^°Hz 


2.1 X 10"^ 
1.9 X 10"^^ 


9.8 X lO^Hz 
3.1 X lO^Hz 


6.0 X 10" 

1.1 X 10" 


-7 

-16 


3.7 X 10"^ < e < 3.8 


X 10"^ 


lO^^GeV 


1.4 X lO'^^Hz 
1.4 X 10^°Hz 


2.1 X lO"'' 
1.9 X 10"^^ 


31Hz 
3.1 X lO^Hz 


2.0 X 10" 
8.2 X 10" 


-7 

-23 


3.7 X 10-^^ < e < 3.^ 


^ X 10"^ 


10^°GeV 


1.4 X lO'^Hz 
1.4 X 10^°Hz 


2.1 X 10-^ 
1.9 X 10-2^ 


3.1 X IQ-^Hz 
3.1 X lO^Hz 


9.5 X 10" 
5.7 X 10" 


-8 
-27 


3.7 X 10-^^ < e < 3.^ 


^ X 10"^ 


lOOGeV 


1.4 X lO'^Hz 
1.4 X 10^°Hz 


2.1 X 10"^ 
1.9 X 10"^^ 


3.1 X 10"'^Hz 
3.1 X 10"^Hz 


5.0 X 10" 
1.3 X 10" 


-9 
-43 


3.7 X 10"^^ < e < 3.8 


X 10-^^ 
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the forecast in Fig. [t] can be a factor of 0(10^^). In order to improve the accuracy of the estimation, we have to run 
simulations with a much larger dynamical range and with many realizations. 




Ti[GeV] 



FIG. 8: The contours of the frequency and the amplitude of the edge in the parameter space of 77 and e (we fixed A = 1.0 and 
= 100). The gray region corresponds t o th e case that the energy density of domain walls dominates the total energy density 
of the universe [the left hand side of Eq. ([28|]. The red region corresponds to the case that e is too large to form domain walls 
with scaling phase [the right hand side of Eq. (28)]. The dotted green line shows the frequency /edge = O.OlHz, IHz, and lOOHz. 
We also show the line with the amplitudes of the edge, Q^^h'^ = 10~^° (orange), 10~^^ (light blue), 10~^^ (pink), 10~^° (blue), 
and 10~^ (light green). 



VI. CONCLUSION 



In this paper, we have studied gravitational waves from domain walls with approximate Z2 symmetry. We have 
performed three dimensional lattice simulations and followed the process of creation, evolution and collapse of domain 
walls. We have confirmed that stable domain walls evolve to maintain the scaling solution, and biased domain walls 
annihilate in the time scale estimated by Eq. (11 ). In numerical simulations we have solved the field equation directly 
rather than using the standard PRS method. It is found that the spectrum of gravitational waves extends over broad 
frequency band, contrary to the naive expectation that gravitational waves have a peak at the frequency corresponding 
to the horizon scale as conjectured in [17, 27 . The results of numerical simulations have revealed that the spectrum 
has the edge corresponding to the horizon scale at the time when walls collapse and the peak corresponding to the 
small scale on the wall typically given by the wall width. In the intermediate range between edge and peak, the 
gravitational wave amplitude increases moderately with frequency, Vt^.j^h? ^ f^-^^. 



By combining analytic estimations and numerical results, we have obtained the formulae (|33[)-(36) and (38) for 
the frequency and the amplitude of gravitational wave spectra. We have extrapolated the results to more generic 
parameter space and showed that the edge of the gravitational wave spectrum from domain walls with energy scale 
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rj = is observable in future gravitational wave experiments such as ET, BBO, and DECIGO. Therefore, 

we expect that the gravitational waves from the collapse of domain walls can be a new observational window to probe 
the models of high energy physics which have (approximate) discrete symmetry. This observation may enable us to 
search the unknown physics through the estimation of the parameters such as bias and provide rich information about 
the theory beyond the standard model of particle physics. 
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Appendix A: Numerical Calculations 



In this appendix we describe the setup for our numerical computations. 



1. Formulation 



We adopt the formulation similar to that of Yamaguchi et al. [49j|50] for the simulation of global strings. The time, 
the Hubble radius, and the temperature are related by the Friedmann equation in the radiation dominated universe 

'=^ = 4-^' (Al) 

where ^ is the constant defined by ^ = (45M|,/167r2^*)^/2 3.68 x 10^^(^*/100)"^/^GeV. We choose the initial time 
ti so that the temperature dX t = ti is twice of the critical temperature of phase transition, Ti = 2Tc = 4/^. We 
normalize the dimensionful quantities in the unit of U. For example, t ^ t/U^ x ^ x/U^ (j) ^ (j)ti, etc. It is useful to 
define the dimensionless quantity ( = ^/rj. By using equation of motion for the scalar field renormalized by ti can 
be written as 

We choose the value of ( as 3.5 which corresponds to = 1.05 x lO^^GeV. 

The (comoving) size of the simulation box is chosen to be b in the unit of t^, and the lattice spacing is Sx = b/N 
where N is the number of grid points. Here we take N = 256. The ratio between the Hubble horizon scale and the 
physical lattice spacing is 

^ = (A4) 

and the ratio between the wall width and the physical lattice spacing is 

-m)-'/\ (A5) 



At the final time t = 151ti, these ratios become /dxpi^ys — 252 < A^, S^/Sxphys — 3.81 for b = 25, and 
i^~^/fephys — 126 < A/", S^/Sxpiiys — 1.90 for b = 50. Therefore, the resolution of the simulations is sufficiently good 
even at the final time. 

We put the periodic boundary condition in the configuration of the scalar field. We solve the time evolution by 
using the fourth order Runge-Kutta method. 



17 



2. Initial Conditions 

As the initial conditions, we assume that the quantum field (j) is in thermal equilibrium with temperature Ti = 
The field and its derivative satisfy the equal-time correlation relations 

(/3|</.(x)0(y)|/3)equal-time = / ^ + (^6) 

— ^[l + 2nky'^<--^\ (A7) 



where Ek = -v/fc^ + and is the occupation number of the Bose-Einstein distributions, ni- = {e^*'/'^^ — 1) ^. The 
mass of the scalar field is given by 



2 d'V 



= \{Tf-Ar,^) = i\r,\ (A8) 

0=0 ^ 



The first term in the square bracket of Eqs. ( A6) and ( A7) corresponds to the vacuum fiuctuations which contribute as 
a divergent term when we perform the integral of k. Hence we subtract this term and use the renormalized correlation 
functions 

(/?|0(x)</.(t/)|/?)re„. = j (0^5- (A9) 

/d^k 
J^E^m. (AlO) 

In the momentum space, these correlation functions can be written as 

{mk)^{k')\/3U^. = ^(27r)35(3)(k + k'), (All) 

{/3|0(k)0(k')|/3)re„. = ^fenfc(27r)3,5(3)(k + k'), (A12) 

where 0(k) is the Fourier transform of (/)(x). Since 0(k) and 0(k) are uncorrelated in the momentum space, we 
generate (^(k) and (^(k) in the momentum space randomly following the Gaussian distribution with 

= ^V, {|^(k)p) = n,E,V, (A13) 

{4>{k)) = iUk)) = 0. (A14) 

Then we transform them into the configuration space and obtain the initial field configurations (/>(x) and (/)(x). Here 
we used (27r)^5*^^^(0) ~ V, where V = b^t^ is the comoving volume of the simulation box. 

3. Calculation of the Area Density 

For calculation of the area density of domain walls, we use the algorithm introduced by Press, Ryden, and 
Spergel [28 . Let us call the neighboring grid points that differ by one in x, z location as the "link". Define 
the quantity S± which takes the value 1 if ^ has different signs at the two ends of a link and the value if has the 
same sign at the two ends of a link. We evaluate the area density A/V hy summing up 6± for each of the links with 
the weighting factor 

A/V = CY^ M II^L U | . (A15) 

where etc. is a derivative of ^(x) with respect to x, and C is a normalization coefficient. We choose C such that 
A/V = 1 is satisfied when all the links have the value 6± = 1. 
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